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We present an algorithm for the retrieval of cloud droplet size distribution parameters (effective radius and 
variance) from the Research Scanning Polarimeter (RSP) measurements. The RSP is an airborne prototype for 
the Aerosol Polarimetery Sensor (APS), which was on-board of the NASA Glory satellite. This instrument 
measures both polarized and total reflectance in 9 spectral channels with central wavelengths ranging 
from 410 to 2260 nm. The cloud droplet size retrievals use the polarized reflectance in the scattering angle 
range between 135° and 165°, where they exhibit the sharply defined structure known as the rain- or 
cloud-bow. The shape of the rainbow is determined mainly by the single scattering properties of cloud par- 
ticles. This significantly simplifies both forward modeling and inversions, while also substantially reducing 
uncertainties caused by the aerosol loading and possible presence of undetected clouds nearby. In this 
study we present the accuracy evaluation of our algorithm based on the results of sensitivity tests performed 
using realistic simulated cloud radiation fields. 

© 2012 Elsevier Inc. All rights reserved. 


1. Introduction 

Understanding of cloud optical and microphysical properties and 
their influences on the interaction of clouds with solar radiation is 
essential for accurate climate projections. Cloud optical properties 
are determined by their liquid water content and droplet size distri- 
bution, which vary within clouds. Accurate and robust remote sensing 
estimates of droplet sizes for different cloud types, especially for bro- 
ken clouds, are crucial for studies of indirect aerosol effects, which are 
focused on relationships between aerosols and cloud properties (such 
as cloud cover, height and particle size). 

Satellite remote sensing methods in the solar spectral domain can 
use both total and polarized reflectance measurements for cloud drop- 
let size retrievals. However, those currently operational are based on 
the multispectral measurements in absorbing and non-absorbing 
bands by e. g., the Moderate Resolution Imaging Spectroradiometer 
(MODIS), and do not include polarization (Nakajima & King, 1990; 
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Platnick et al., 2003). These retrievals have certain limitations owing 
to the influence of the 3D structure of clouds (especially broken cumu- 
lus and stratocumulus) on the reflected radiation, which is not 
accounted for in the ID radiative transfer models used in the retrieval 
algorithms (Marshak et al., 2006; Zinner et al., 2010). In particular, 
Marshak et al. (2006) found that ignoring sub-pixel cloud variability 
produces a negative bias in the retrieved effective droplet radius, 
while ignoring cloud inhomogeneity at scales larger than a pixel 
scale, on the contrary, leads to overestimation of the average droplet 
size. The 3D shadowing and illumination effects also introduce depen- 
dence of the retrieval results on solar and viewing geometries. 
Vant-Hull et al. (2007) found that when 3D clouds are viewed near 
the backscatter geometry, the retrievals based on a plane-parallel 
model underestimate the effective radius, while the reverse is true 
when the satellite is far from the backscatter position. The satellite 
data survey by Girolamo et al. (2010) found that the retrieved cloud 
reflectance is indistinguishable from plane-parallel clouds in only 
24% of cases, mostly limited to regions dominated by stratiform clouds 
at solar zenith angles less than 60° (for other regions or solar angles 
this frequency drops sharply to as low as a few percent). Nakajima 
et al. (2010) found that the cloud droplet sizes retrieved from 
MODIS data would be strongly influenced by the vertical inhomogene- 
ity of droplet size, hypothesizing the existence of small cloud droplets 
at cloud top as the reason for small droplet size retrievals using the 
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3.7-pm channel and large drizzle drops influencing 2.1- and 1.6-pm 
channels. The sensitivity of MODIS particle size retrievals to drizzle 
formation was not, however, corroborated by Zinner et al. (2010). 
The gaseous (water vapor) absorption inside the cloud also creates 
uncertainties (positive bias) in retrievals of water cloud effective radi- 
us from total reflectance measurements (Platnick& Valero, 1995) and 
a similar effect can result from the presence of an absorbing aerosol 
layer above the cloud (Coddington et al., 2010; Haywood et al., 2004). 

Retrieval of cloud droplet size from polarized observations of the 
reflected light in the rainbow region (at scattering angles between 
135° and 165°) allows many of the uncertainties associated with 3D 
effects and unknown aerosol loadings to be avoided. This advantage 
results from the shape of the rainbow being dominated by single 
scattering of light by cloud particles, even though its amplitude can be 
affected by aerosols and the geometric structure of clouds. The domi- 
nance of single scattering as a determinant of the structure of the rain- 
bow simplifies both forward modeling and inversions. The retrievals are 
just as accurate over land, or ocean (no surface albedo issues), and are 
valid independent of the optical depth down to unity (i.e. they work 
for common low water path clouds). 

Another potential use of this technique is in combination with lidar 
measurements of cloud-top extinction. The polarimetric retrievals 
of effective radius and variance provide accurate estimates of the 
extinction cross-section of cloud droplets that, when divided into the 
extinction coefficient, allow the number concentration of cloud droplets 
to be determined. This quantity, unlike the optical characteristics, is 
directly related to aerosol effects on clouds (Brenguier et al., 2003). 
Hu et al. (2007) introduced a method for derivation of droplet number 
concentration from CALIPSO satellite lidar data combined with MODIS 
droplet size retrievals. Their approach used an empirical relationship 
between depolarization and cloud extinction coefficient (obtained 
using Monte Carlo simulations) and derived the extinction cross- 
section from the MODIS 3.7-pm channel size retrieval. However, since 
the MODIS dataset does not provide effective variance of the droplet 
size distribution there are non-negligible uncertainties in the estimate 
of extinction cross-section that affect the droplet number concentra- 
tion. This obstacle would not affect a study where polarimetric 
retrievals of droplet size distribution are used. 

The polarized rainbow technique has previously been used to 
retrieve cloud droplet effective radii from the Polarization and Direc- 
tionality of the Earth's Reflectances (POLDER, Deschamps et al., 1994) 
measurements (Breon & Doutriaux-Boucher, 2005; Breon & Goloub, 
1998). The inter-comparison study by Breon and Doutriaux-Boucher 
(2005) found high correlation between POLDER and MODIS products 
over the oceans (with MODIS biased high by about 2 pm). This corre- 
lation, however, breaks down over most of the continental and 
polluted oceanic areas, where POLDER retrieves smaller droplets 
(less than 7 pm effective radius). 

While RSP allows both total and polarized reflectance methods to 
be used simultaneously, in this study we focus only on the polarized 
rainbow technique. Here we present the results of tests performed 
on the output from ID and 3D radiative transfer models that are 
used to simulate the polarized reflectance of both simplified and real- 
istic cloud scenes. Our retrievals from actual RSP measurements 
performed during two field campaigns (to be published elsewhere) 
appear to be consistent with the correlative in situ observations 
within 50-100 m of the cloud top. 

2. RSP measurements 

The RSP (Cairns et al., 1999) is an airborne prototype for the satel- 
lite Aerosol Polarimetery Sensor (APS), which was built as part of the 
NASA Glory Project (Mishchenko, 2006; Mishchenko et al., 2007). 
This instrument measures /, tQ, and U components of the Stokes 
vector (cf. Hansen & Travis, 1974; Mishchenko et al., 2006) in 9 spec- 
tral channels with center wavelengths of 410, 470, 555, 670, 865, 960, 


1590, 1880 and 2260 nm. The total and polarized reflectances (R and 
R p respectively) are then derived from these Stokes parameters as 


and 


r p = 


nQ 


( 2 ) 


where I 0 is the extraterrestrial solar irradiance, and jM, is the cosine of 
the solar zenith angle (SZA). The Stokes parameter Q in Eq. (2) is 
defined with respect to the scattering plane containing both solar 
and view directions (parameter U in this plane is negligibly small). 
In our notation the polarized reflectance is positive when the polari- 
zation direction is parallel to the scattering plane and negative 
when it is orthogonal to that plane. (Note, that this notation is the 
opposite to that adopted by Breon and Goloub (1998) and Breon 
and Doutriaux-Boucher (2005). 

The RSP scans its 14 mrad field of view over ±60° from nadir 
(starting at forward direction) and making samples at 0.8° intervals, 
with additional calibration measurements being made in other parts 
of the scan. Thus, the observational data obtained in each scan 
consists of around 150 instantaneous Earth viewing measurements. 
Usually the RSP is oriented to scan along the aircraft track. In this con- 
figuration the data from the actual RSP scans can be aggregated into 
“virtual” scans consisting of the reflectance at the full range of view- 
ing angles at a single point on the ground, or at the cloud top 
(Fig. 1 ). This aggregation can be done using either the aircraft attitude 
data (altitude, speed, pitch and crab angles), or the RSP measure- 
ments of brightness contrast between small clouds and the surface. 

In this study we focus on liquid water clouds with cloud tops 
warmer than 0 °C. For colder cloud tops, the cloud thermodynamic 
phase (ice, water or mixed-phase clouds) can be determined by ana- 
lyzing the magnitude of the main rainbow feature (Goloub et al., 
2000; van Diedenhoven et al., 2012b). Besides liquid water cloud 
properties RSP measurements can be used for accurate retrievals of 
ice cloud parameters (Ottaviani et al., 2012; van Diedenhoven et al., 
2012a), aerosol optical depth, size, and refractive index (Cairns et 
al., 2009; Knobelspiesse et al„ 2011a, 2011b; Waquet et al., 2009), 
as well as for characterization of the ground surface (Knobelspiesse 
et al., 2008). Over ocean, in addition to aerosol retrievals, chlorophyll 
concentrations are simultaneously estimated from the atmospheri- 
cally corrected ocean color measured by the RSP (Chowdhary et al., 
2005, 2006, 2012; Gilerson et al., 2006). 


3. Retrieval algorithm 

For cloud droplet size retrievals we use the scattering angle 
dependences of the polarized reflectance. In this study we mostly 
use the 865 nm spectral channel, however, other RSP channels can 
also be used (except for 960 and 1880 nm channels affected by 
water vapor and/or C0 2 absorption). Our technique is focused on 
the sharply defined structure (rainbow) in the polarized reflectance 
of clouds within the scattering angle range between 137° and 165° 
(Fig. 2). The dependence of the rainbow signature on the cloud drop- 
let effective radius has the form of a dilation of the curve along the 
scattering angle axis, while increase of the effective variance results 
in smoothing of the curve making the extrema less pronounced 
(until their eventual disappearance). 

To analyze the shape of the polarized rainbow, we first perform a 
rotation from the measurement coordinate frame to that of the scat- 
tering plane (cf. Hansen & Travis, 1974). Since the rainbow shape is 
determined mainly by the single scattering, this procedure results 
in an increase in magnitude of the Stokes parameter Q, while the 
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Fig. 1. Schematic illustration of the RSP data aggregation process: the data from the actual scans is aggregated into "virtual'' scans, each consisting of all reflectances (at different 
viewing angles) from a single point at the cloud top. 


parameter U becomes an order of magnitude smaller than Q (U would 
become 0 in the case of complete absence of multiple scattering). 
Then, we fit the polarized reflectance computed according to Eq. (2) 
by functions from the following parametric family: 

Kp(y) = :A'C(T + 6; r ef f , v eff ) + Bcos 2 y + C. (3) 

Here y is the scattering angle, r eS and v e fr are respectively the 
effective radius and variance of the cloud droplet size distribution 
(cf. Hansen & Travis, 1974), which is assumed to have the form of a 
gamma distribution. The phase matrix elements pj“ le) are computed 
using Mie theory for a grid of effective radius and variance values, 
while A, B, and C are empirical fitting parameters accounting for contri- 
butions to the polarized reflectance from everything beyond single 
scattering by cloud droplets, such as multiple scattering, Rayleigh scat- 
tering, aerosol extinction, ground surface reflectance for thin clouds, as 
well as effects caused by rotation to the scattering plane. Furthermore, 
the effect of overlying cirrus clouds is expected to be largely mitigated 



by incorporating these empirical fitting parameters in the procedure 
(B in particular), since measurements indicate that natural ice clouds 
generally have featureless polarized reflectances, which are near- 
linear as a functions of scattering angle (C.-Labonnote et al., 2000; 
Knap et al., 2005; van Diedenhoven et al., 2012b). Small fitting parame- 
ter 6 is supposed to compensate for angular errors (shifts) in measured 
and simulated reflectances. The cos 2 y term is included to capture any 
Rayleigh scattering contributions to the observations. Our tests showed 
that this parameterization provides a slightly better fit to simulated data 
than the parameterization of Breon and Goloub (1998), which uses a 
term linear in y. Note, however, that the difference between the two 
parameterizations is small since in the rainbow region cos 2 y can be 
closely approximated by a linear function of y (plus a constant). 

We use a pre-calculated look-up table (LUT) of values with 
0.2° resolution in the scattering angle. The grid values of r eff used in 
this LUT range from 5 to 20 pm with 0.5 pm increments, while the 
grid for v eff starts with the values 0.01, 0.03, 0.05, and then becomes 
regular with 0.025 increments and maximum value of 0.35. Our 
fitting technique consists of 2 steps: first, we count minima and 



Fig. 2. Sensitivity of polarized reflectance in the rainbow angular range to effective radius (left) and variance (right) of the cloud droplet size distribution. The data are simulated for 
865 nm wavelength in single scattering approximation. 
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maxima in the observed rainbow signature and match these numbers 
(give or take one) with those from the corresponding LUT; after that 
we directly look for the best fit (according to Eq. 3) among the plau- 
sible subset of forward modeled curves selected on the first step. 
After the best fit cloud size model is determined we perform a refine- 
ment procedure in its neighborhood on a 10 times denser grid in 
effective radius and variance. The first step can be skipped for analysis 
of simulated data, unless noise and/or measurement artifacts are 
deliberately introduced into the simulations. 

4. Tests on ID plane-parallel simulations 

To justify our (semi-empirical) method, we need to demonstrate on 
simulated data that the parameterization (3) adequately separates the 
cloud single scattering contribution from other factors providing high 
accuracy retrievals. We use a vector multiple scattering radiative trans- 
fer (RT) model to simulate RSP measurements of polarized reflectances 
emerging from a cloud field with known microphysical parameters. 
Then we derive cloud droplet size from these data using our retrieval 
algorithm, and compare the retrieved values of the droplet effective 
radius and variance with those assumed in the simulations. 

The modified vector doubling/adding code (Cairns et al., 1997) 
used for forward modeling of reflectances for plane-parallel (ID) 
atmospheres is essentially the same as that presented by Hansen 
and Travis (1974) with modifications based on the work of de Haan 
et al. (1987). The simplicity of this RT model allowed us to perform 
several numerical experiments in order to study the sensitivity of 
the retrieval accuracy to various factors, such as the scene geometry, 
presence of a thick aerosol layer above the cloud, and composition of 
the cloud itself. We do not investigate the influence of the ground sur- 
face, since its polarized reflectance is generally smooth as a function 
of scattering angle and should be eliminated by the regression (3). 
All simulations described below were made for the 865 nm wave- 
length, unless specified otherwise. 

We start with the most straight-forward numerical experiment, 
when the RSP data is simulated for the measurements in the principal 
plane (containing the Sun direction and the normal to the surface). In 
this case the Sun is directly ahead or behind the aircraft (so that the 
relative solar azimuth is either 0° or 180°). The RSP measurements 
in this case cover the maximal scattering angle range for the given 
SZA. No transformation of the Stokes vector is needed, since in this 
situation the scattering plane corresponding to each viewing angle 
coincides with the principal plane. Thus, no errors are introduced in 
association with such a transformation. We performed this experi- 
ment on an extensive grid of realistic cloud size parameters, with 
r e ff values of 5, 7.5, 10, 12.5, 15, and 17.5 pm, and v e ff values of 0.01, 
0.05, 0.1, and 0.2. In all cases SZA was 60°, and the cloud optical 
depth (COD) was five, which is effectively infinite for the generation 
of polarization. We allowed the angular shift fitting parameter S to 
vary within ±0.2° (which is the resolution of both the forward 
model and the LUT used for inversions) with 0.01° increments. The 
fitting errors are computed as the standard deviations of the differ- 
ences between the simulated and fitting values on the angular grid. 
Examples of the simulated and best fit (according to Eq. (3)) polar- 
ized reflectance are shown in Fig. 3 by respectively red and green 
curves for r eff =7.5, 17.5 pm, and v eff =0.01, 0.2. The blue curves 
show the reflectance computed assuming exactly the same size distri- 
bution parameters as those used in the RT model. They are slightly 
different from the best-fit curves, partly due to the small angular un- 
certainties (best seen in the top right panel). The comparison between 
the retrieved effective radius and variance values with those assumed in 
the simulations is shown in Fig. 4 as separate scatter plots for r eff and v e ir 
and also as the difference vectors in (r eff , v e ff)-plane. There the dia- 
monds represent the assumed models, and the triangles are the 
best-fit retrievals. We see from Fig. 4 (top) that the mean accuracy of 
the r eff retrievals is better than 0.1 pm, while its standard deviation 


has maximal value of 0.21 pm at v eff =0.2 (due to the single worst 
case shown in Fig. 3 (top right)). While the error in the effective radius 
retrievals is rather random, the derived values of v eff show systematic 
biases (Fig. 4 (middle)), which can be approximated by r e frdependent 
linear functions of the assumed v eff (they are listed at the bottom of 
the plot). We see that these biases range from 6% to 27% and generally 
decrease with r eff , while the largest absolute differences occur at large 
v e fr, especially 0.2. This may be explained by the lack of “structure" 
(number and amplitude of oscillations) in the angular dependence of 
polarized reflectance for both small sizes and large variances, as can 
be seen in Fig. 2. The “smoothing” effect of multiple scattering on the 
polarized reflectance, is then interpreted by the retrieval algorithm as 
a bias toward larger v e ff, and (to a much lesser extent) to smaller droplet 
size. 

Similar small biases can be observed in droplet sizes retrieved 
from different spectral channels. Fig. 5 shows polarized reflectances 
and size retrievals for the same models as presented in the top panels 
of Fig. 3 (r e ff=7.5 pm, v e ff=0.01, 0.2), but for a shorter (410 nm) and 
a longer (1588 nm) wavelength. The increasing lack of “structure” in 
these curves causes some underestimation (by 0. 1-0.6 pm) of the 
droplet r eff as wavelength increases, and the larger retrieval errors 
are observed for large v e ff. While these biases are relatively small for 
simulated data, the presence of measurement noise or artifacts is like- 
ly to more strongly affect the retrievals when the rainbow signature 
is weak. Thus, we do not recommend using the longest RSP wave- 
lengths (1588 and 2260 nm) alone. Our tests show that polarized 
reflectance in a single channel is sufficient for accurate determination 
of droplet size distribution parameters. However, simultaneous 
retrievals using multiple spectral channels are beneficial for providing 
a retrieval consistency check, as well as a quality assessment of the 
measurements (particularly for airborne measurements where 
pointing knowledge is unlikely to be better than ±1°). 

5. Effects of rotation to the scattering plane 

Before describing the effect of the coordinate frame rotation on the 
droplet size retrievals, we should note that the range of the observed 
scattering angles also depends on solar/viewing geometry, and this is 
an important factor affecting the retrieval accuracy. Obviously, good 
retrievals cannot be expected, when only a small part of the rainbow 
structure is seen in the data. Fig. 6 shows the dependences of the scat- 
tering angle range on the relative azimuth between the Sun and aircraft 
directions for three different solar zenith angles (40°, 60°, and 80°), 
assuming that the viewing angle varies within ± 60° from nadir. The 
relative azimuth intervals are classified in these plots according to the 
extent (“Full”, “Partial”, or “None”) to which the rainbow region (here 
140°-165°) is covered by the RSP measurements made under the 
corresponding viewing geometries. These results are summarized in 
the diagram in Fig. 6 (bottom right). We see from this diagram that in 
the limit case, when the Sun is directly at zenith, the full rainbow 
range is observed independently of the relative Sun/aircraft azimuth. 
For non-zero SZA the best rainbow coverage is achieved in the principle 
plane (zero relative azimuth), then, the coverage decreases with the 
azimuth as it increases to 90° (values larger than 90° are not shown 
due to symmetry). The coverage also gets worse with increase in SZA, 
so that for SZA>75° full-range rainbow observation cannot be made 
even in the principal plane while at SZA >40° there is a relative azimuth 
range, where the rainbow cannot be seen at all. 

As we already mentioned, the polarized reflectance measured in 
the solar principal plane is determined by the Stokes parameter Q 
according to Eq. (2), while the parameter U is zero as a result of sym- 
metry. In order to perform the droplet size retrievals in an arbitrary 
viewing geometry the Stokes vector needs to be transformed from 
the measurement coordinate frame into that of the scattering plane. 
Note, that this plane is defined as the one containing both solar and 
viewing directions, thus, in a generic case it is different for different 
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Fig. 3. Retrievals of effective radius and variance from the polarized reflectances simulated using plane-parallel adding/doubling vector code for r e ff=7.5 and 17.5 pm, and 
v e(f =0.01 and 0.2. The simulated (using multiple scattering) reflectances are shown in red; the single-scattering fits assuming the same r eff and v eff that were used in the radiative 
transfer model are in blue; while the best fits (assuming an angular shifts) are in shown in green. In all cases COD is 5. 


viewing angles within the same RSP scan. In the case of single scatter- 
ing, the transformation of the Stokes vector corresponding to the 
rotation of the reference frame is linear (cf. Hansen & Travis, 1974) 
and affects only its Q and U components (U=0 after this rotation). 
While, strictly speaking, such a simple transformation is not valid in 
the case of multiple scattering in clouds it can still be used in the rain- 
bow region, where the angular structure of polarized reflectance is 
dominated by single scattering. In this case the transformed U is 
orders of magnitude smaller than the transformed Q which has simi- 
lar structure in the rainbow to that of single scattering. The deviation 
in the parameter Q from that of single scattering, resulting from 
the reference frame rotation under multiple scattering conditions, is 
expected to be a smooth function of scattering angle, which is 
absorbed by the last two terms in the parameterization (3). This 
assumption is confirmed by the excellent fits of the rotated polarized 
reflectances using Eq. (3) in the numerical tests described below. 

We applied our retrieval algorithm to RSP measurements simulat- 
ed for 3 solar zenith angles (20°, 40°, and 60°) and the corresponding 
relative azimuths (48°, 22°, and 16°) chosen to be the most distant 
from 0°, while still allowing for full rainbow range observation (see 
Fig. 6). This choice was made to separate the errors introduced by 
the frame rotation from the loss of accuracy due to reduced angular 
range. For comparison we also generated and analyzed data for the 
same SZAs, but with the measurements made in the principal plane. 
The cloud parameters assumed in simulations were the same in all 
cases: r eff =17.5 pm and v e fr=0.1,and COD = 5. The polarized reflec- 
tances obtained after the transformation into scattering plane frame 
are shown in Fig. 7. We see that the “rainbow signal” is weaker at 


smaller SZA, introducing larger uncertainties in the retrievals, espe- 
cially for SZA = 20° (which, however are still quite good: only 
0.5 pm underestimation of r e fr, and 10% overestimation of v e fr). The 
“Lf-reflectances” computed using the analog of Eq. (2), but with Q re- 
placed by the residual value of U in the scattering plane are shown in 
Fig. 7 (bottom right). One can see that these residues are more than 2 
orders of magnitude smaller than the corresponding polarized reflec- 
tances (while larger than the corresponding quantities in the princi- 
ple plane simulations, which are essentially 0 (smaller than 1CT 9 )). 
The differences between retrievals from the “rotated" datasets and 
from their principle plane counterparts (r e ff=17.2 pm, v e fr=0.110 
for SZA = 20°, r eff =17.5pm, v eff =0.107 for SZA = 40°, and r eff = 
17.5 pm, v e ff= 0.102 for SZA = 60°) are as little as 0.1-0.2 pm in radi- 
us and a few percentage points in variance. This allows us to conclude 
that the error introduced by the transformation of the Stokes vector 
to the scattering plane frame and omission of U is negligible. 

6. Effects of angular shifts 

Angular errors, if unaccounted for, can significantly affect the 
accuracy of satellite, and especially airborne measurements. In partic- 
ular, inadequately recorded aircraft pitch, roll and yaw may cause an 
error in the computed scattering angle. In realistic situation this error 
is a non-linear function of the aircraft attitude uncertainties, which 
can be large. However, if these uncertainties are small (less than a de- 
gree), they result in an effective shift in scattering angle, which is easy 
to model. To see the effect of such a shift on the droplet size retrievals, 
we took the polarized reflectance simulated for clouds with r efr = 7.5 
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Fig. 4. Comparison of retrievals of effective radius (top) and variance (middle) with 
the actual values assumed in plane-parallel vector simulations (in all cases COD = 5). 
Bottom: the same comparison in ( v e ff)-plane. Diamonds depict the actual values 
assumed in simulation, triangles — best fit retrievals. 


and 17.5 pm and displaced the actual scattering angle records by 
small angle Ay ranging from -0.5° to 0.5°. In all considered cases 
v e ff = 0.1, COD = 5, SZA is 60°, and the relative solar-aircraft azimuth 
is 0°. The results of application of our algorithm (without optimiza- 
tion for angular shift parameter 6) to this dataset are presented in 
Fig. 8 for the case of r ef f= 17.5 (plots in the case of r e ff = 7.5 are sim- 
ilar). A strong linear dependence of the retrieved r eff on Ay is seen in 
the top right panel of Fig. 8. This seems reasonable when one con- 
siders that, while the dependence of the rainbow shape on the drop- 
let effective radius can be characterized as a dilation with fixed point 
at 135° in scattering angle space (as r eff decreases), small shifts can 
also be interpreted as dilations. To see this, we can apply shift Ay 
to both ends of the scattering angle interval [y 0 ,yi]- At each end 
the shift can be interpreted as a multiplication by a factor c,- = 
1 + (Ay/y,), i = 0.1. The factors c 0 and Ci are identical in the case of 
true dilation, while in our case they are very close as Ay<Ky 0 <yi. 

The accuracy of the numerical simulations used for these tests was 
0.2°, so it is not surprising to see a discrepancy of that order between 
simulated multiple- and single-scattering polarized reflectances. 
In the case of r eff = 7.5 pm the size assumed in the model can be re- 
trieved when applying a -0.25° shift, while the effective radius 
value retrieved without a shift was 7.2 pm. Given the linearity of 
the dependence of the retrieved r eff on the implied Ay, this means 
that the retrieval error in this case is 1.2 pm or 16% per degree of 
the shift. In the case of r eff = 17.5 pm (Fig. 8) a smaller —0.15° shift 
was necessary to recover the assumed droplet size, while r e «= 
16.8 pm was retrieved with Ay = 0. This corresponds to 4.7 pm or 
27% error rate per degree of the shift. The error associated with angu- 
lar uncertainty is more significant for larger particles, since the oscil- 
lating structure of the rainbow generated by them has shorter angular 
period than when the droplets are smaller (see Fig. 2). Thus, the same 
angular shift appears to be larger relative to this period for larger 
droplets, and, therefore, has stronger effect on retrievals. The effect 
of the shift on the retrieved v eff is weaker, and practically negligible 
in the case of r eff = 17.5 pm (Fig. 8 (middle right)). 

The bottom right panel of Fig. 8 demonstrates that the uncer- 
tainties in scattering angle shift can be effectively resolved if the mag- 
nitude of the shift is taken as a free parameter in the fitting procedure, 
since the best fits to the simulated reflectances occur with the same 
Ay, for which the retrieved values of r eff coincide with those assumed 
in the model. This feature is implemented in our algorithm by intro- 
duction of the shift parameter 5 in Eq. (3), which is allowed to vary 
within ±0.2° for simulated data tests. In the case of real airborne 
data analysis, when the measurements can be affected by uncer- 
tainties in aircraft attitude, we first correct the pitch angle, so that 
the first (and the deepest) minimum in polarized reflectance is 
between 140° and 145°, and then apply fitting procedure, which 
includes optimization with respect to smaller angular errors. 

7. Aerosol layer over cloud 

We performed a number of simulations of RSP measurements over 
a mixed atmosphere consisting of a cloud and an aerosol layer on top 
of it. This allowed us to study the effect of aerosols on cloud droplet 
size retrievals, and also to evaluate the RSP algorithm's sensitivity to 
a rainbow signal originated deep in the scattering medium. All simu- 
lations were performed for the same viewing geometry (principal 
plane, SZA = 60°) and cloud parameters (r efr = 17.5 pm, v eff =0.1, 
and COD = 5). The aerosol layer was modeled by small water particles 
(n r = 1.3275 (same as in cloud), r eff =0.15 pm and v e fr=0.1), and had 
varying optical depth (AOD = 0, 0.1, 0.25, 0.5, 0.75, 1, and 2). The 
polarized reflectance of the aerosol alone increases monotonically 
with the scattering angle between 137° and 165°, and does not exhib- 
it any rainbow features. Addition of the aerosol layer on top of the 
cloud leads to a decrease in the “rainbow signal" magnitude, that 
affects the droplet size retrieval accuracy. The dependence of the 
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Fig. 5. Same as Fig. 3 (top) but for 410 nm (top) and 1588 nm (bottom) wavelengths. 


retrieved droplet effective radius and variance on the thickness of the 
aerosol layer is weak, and the retrieval errors are similar to the 
cloud-only case: less than 0.5 pm in r eff and 0.04 v e ff. The largest er- 
rors occur at the largest AOD value (0.75), for which retrievals were 
still possible despite an increasing lack of “structure” in the dataset. 

Fig. 9 illustrates a simple test characterizing the effect of the aero- 
sol layer on the rainbow signal magnitude, based on the simulations 
described above. We took all polarized reflectances simulated for 
different AOD (top left panel in Fig. 9) and represented them in the 
form similar to Eq. (3) 

Kp(y; Ta) = AR p (y; = 0) + By + C, (4) 

where r a is AOD, and A'=A'(r a ) is the magnitude of the rainbow 
signal relative to the case with no aerosol presence (r a = 0). Unlike 
Eq. (3), this form does not imply droplet size retrievals or modeling 
using Mie theory (we use prime accents for the coefficients to empha- 
size this difference). The model fit (4) is good and the residues 
^(yifa) — /A'RpCyjO) can be very closely fitted by linear functions of 
the scattering angle y (Fig. 9(top right)). These residuals can be 
attributed to the polarized reflectance of the aerosol layer. The depen- 
dence of the signal magnitude A' on AOD appears to be exponential 
(see Fig. 9 (bottom left)), and agrees well with the Beer's law attenu- 
ation model: 



Here Ps and are cosines of respectively the solar zenith angle d s 
and the viewing angle 0 V . Knowing the solar zenith angle (0 S = 6O°) 
used in our simulations, we can estimate (see Fig. 9 (bottom right)) 


the value of the average viewing angle (0 V = 29°) from Eq. (5), and 
using that 

y = 180 —0 S + 0 V , (6) 

find the corresponding scattering angle y=149°. As one would 
expect, this value is in the middle of the rainbow angular range. 
Thus, our numerical experiment shows that the magnitude of the 
rainbow, in polarized reflectance observations, is attenuated by over- 
lying aerosols according to the Beer’s law but its structure is not 
changed by the overlying aerosol scattering. 

8. Contribution of forward scattering to the rainbow 

Our tests described in this section show that the polarized rainbow 
amplitude includes contributions of both single and forward-directed 
multiple scattering. The latter, does not affect the angular structure of 
the rainbow, and, therefore, droplet size retrievals. However, assess- 
ment of the magnitudes of various contributions will help to under- 
stand the retrievals from stratified clouds described in detail in the 
next section. 

Let us determine the general functional form that the reflectance R 
of a single-layer cloud of optical depth r should take in the single- or 
forward-scattering case. R must satisfy the following property, when 
an additional layer of COD = is attached to the top of the cloud of 
COD = t 2 : 

R(Ti+t 2 ) = R(Ti) + u(t 1 )R(t 2 ). (7) 

Here u(r) = exp(— j3r) is the attenuation factor. The exponent 0 
depends on the scattering geometry and may also be influenced by 
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Fig. 6. Top, and bottom left: dependences of the scattering angle range viewed by RSP (light gray) on the solar azimuth (relative to the aircraft orientation) for three values of the 
solar zenith angle. The angular ranges where rainbow can be observed (fully or partially) are shown by darker shades of gray. Bottom right: diagram showing availability of obser- 
vations in the rainbow angular range (140°-165°) depending on the solar zenith angle and the relative solar azimuth. 


multiple scattering in the forward direction, which depends on the 
droplet size. If we look at a rather narrow scattering range (as in the 
rainbow case) the explicit dependence of the attenuation factor u on 
the scattering angle can be neglected, and we can assume that the 
reflectance R can be factorized into a product of the r-independent 
angular spectrum and an attenuation function that depends on the 
optical depth only. Accepting these assumptions, it is not difficult to 
show that the only functional form of R(r) satisfying the group transfor- 
mation property (efeq: R — tl + f2) is 

R(r) = [l-e _,3T ] (8) 

where R ( ” ) is the reflectance of semi-infinite medium with the same mi- 
crophysical properties. 

To test this formula, we simulated the polarized reflectance for 
clouds with droplet r eff = 5, 7.5, 12.5 and 17.5 pm (in ail cases v eff = 
0.1, SZA = 60°, measurements are made in the principal plane). COD 
in these tests ranged between 0.1 and 0.5 with 0.1 increments. R M 
was taken from simulations with COD = 5, which is essentially infi- 
nite for the generation of polarization, in all cases Eq. (8) was satisfied 
with very high accuracy. This is demonstrated in Fig. 10 for the case of 
r eff = 17.5 pm. The polarized reflectances generated for different COD 
values were fitted by scaled R M using procedure similar to Eq. (4). 
The corresponding scaling factors v(r) are then used to compute the 
attenuation functions u(r) = 1 — v(r), which appeared to have an 
exponential form. The values of attenuation exponents corresponding 
to r rm eff= 5, 7.5, 12.5 and 17.5 pm are respectively 0 = 1.91, 1.89, 1.79, 
and 1.57. All these numbers are smaller than the single-scattering 


airmass 0 = 3.14, indicating that the polarized rainbow amplitude 
includes contribution from forward scattering within the phase func- 
tion diffraction peak. This is also supported by the fact that the 
observed 0 is smaller for larger particles, that have stronger forward 
scattering. At the particle size close to that of aerosol (0.15 pm), 0 is 
back to the single-scattering value (as it was shown in the previous 
section). 

9. Stratified cloud 

The size distribution of droplets in real clouds usually varies with 
the altitude. This raises the question of how deep can RSP “see” into a 
cloud and what its measurements represent regarding the statistics 
of the cloud droplet size profile. The assumption, that the polarized 
rainbow shape is dominated by single scattering, allows us to better 
understand the relationship between the droplet size distribution pa- 
rameters retrieved from polarized reflectance at the cloud top and the 
profiles of these parameters inside the cloud. 

9.1. General properties 

Let us consider a vertically stratified cloud of the optical depth 
Tmax. and parameterize the altitude dependence of droplet size distri- 
bution in it by the optical depth r measured from the top of the cloud 
downwards (0< r<r max ). Suppose that we successfully extracted the 
single-scattering component PS^Hy) of the polarized reflectance 
measured by RSP at the cloud top. We assume that the cloud particles 
contribute to PS-p^y) independently one from another, thus, it is a 



100 


M.D. Alexandrov et al. / Remote Sensing of Environment 125 (2012) 92-1 1 1 



Fig. 7. Top, and bottom left: Polarized reflectances (after transformation to the scattering plane) and the corresponding retrievals from RSP measurements. The data are simulated 
for different Sun/viewing geometries: SZA of 20°, 40°, and 60°, and the maximal relative solar-aircraft azimuth still allowing for full rainbow range observation at given SZA. In all 
cases r e ff= 17.5 pm and v e ff=0.1, and COD = 5. Bottom right: the residual “l/-reflectances” left after transformation to the scattering plane in these simulations. 


superposition of the signals coming from different depths. This allows 
us to introduce the “effective” droplet size distribution n top (r). This 
distribution is derived from the RSP measurements, since it deter- 
mines the polarized reflectance at the cloud top: 

J ^sca (r)P u (y^)n 

top (r)dr 

P ( u ' P) (Y) = * (9) 

J'o'sca(r)n top (r)dr 

o 

where cr sca (r) is the scattering cross-section. The cloud-top distribu- 
tion ritop(r) is normalized by the condition 


and n r (r) is the droplet size distribution in the optical depth interval 
[r,r+ dr]. If the latter distribution is assumed to be normalized by the 
condition 

Jn r (r)dr = l, (13) 

o 

its weight w(t) in n top (r) is determined by two factors: the number of 
particles dJV c (per unit area) in the corresponding optical depth inter- 
val dr, and the Beer's law attenuation factor u(r): 

w(r) = Cu(r) t ^. (14) 


J n to P ( r )dt = 1, (10) 

o 

and it can be represented as a superposition of size distributions at 
different optical depths: 

T max 

n to P ( r ) = J w(r)n T (r)dr, ( 11 ) 

0 

where w(t) is the weighting function normalized as 

T max 

j w(r)dr = 1 . ( 12 ) 

o 


Here C is the normalization constant derived from Eq. (12), and 


u(t) = exp[— fir] = exp 


~ ( — + — ] CLT 

A 


(15) 


where a provides a correction for the effects of multiple scattering in the 
forward direction. The extinction optical depth of a homogeneous layer is 

t = N c (a ext ) = N c (nr 2 Q ext (r )j = N c rr(Qext>(r 2 ) (16) 


with N c being the 2D column particle concentration (i.e. the number of 
particles in a column over a unit area), cr ext and Qe Xt are respectively 
the extinction cross-section and efficiency in the layer, and the averaging 
is over a statistical ensemble of particle sizes specified by the size 
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Fig. 8. Left: Polarized reflectances and retrievals from simulated RSP measurements with artificial scattering angle shifts: Ay= -0.2°(top), 0°(middle), and 0.2°(bottom). In all cases 
r teff = 17.5 pm, v eff =0.1, COD = 5, SZA is 60°, and the relative solar-aircraft azimuth is 0°. Right: Sensitivity of the retrievals from simulated RSP measurements to scattering angle 
shifts ranging from —0.5° to 0.5°. 


distribution. The average extinction efficiency<Qext>is defined by 
Hansen and Travis (1974) as 


(Qext) 


G 


(nr 2 Qex t( r )) 

W) 


(17) 


where G is the average geometric cross-section. It follows from Eq. (16) 
that 


specific droplet size. If particle radius is measured in microns, the 
units of the concentration N c will be pm~ 2 . 

While the distribution n top (r) can have an arbitrary functional 
shape, we can estimate its effective radius 


.-(top) 



(19) 


^ = 1 , (18) 

dT men)r(r- 2 ) T 

wher e<r 2 > T is the second moment of the particle size distribution 
n T {r). Note, that for large cloud particles <Q ext >« 2 regardless of the 


and variance 


„< t0 P> 



- 1 , 


( 20 ) 
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Fig. 9. Summary of the retrievals from the mixed atmosphere simulations. Top left: polarized reflectances vs. scattering angle for a variety of the top layer AOD. Top right: the 
residues of the reflectances from top left panel remaining after subtraction of the (scaled with the fitting factor a) reflectance corresponding to AOD = 0 (solid lines), and linear 
fits to these residues (dashed lines). Bottom left: dependence of the fitting factor a on AOD. Bottom right: inversion of the mean scattering angle from this dependence. 


which are the primary determinants of the cloud's optical properties and a straight-forward calculation shows that 
(cf. Hansen & Travis, 1974). The moments of the distribution ( 1 1 ) are 


(Otop = J r m n top (r)dr = J <r m ) T w(r)dr, 
0 0 


( 21 ) 


..(top) 



where<r m > T are the corresponding moments of the distribution 

Mr). 

Let us assume for simplicity that n r (r) has gamma distribution 
shape (cf. Hansen & Travis, 1974) 


and 



n(r) = ^ (2b ~ m 
{> r[(l— 2b)/b] 


r (l— 3 b)/b g -r/ab 


( 22 ) 

9.2. 2-layer cloud 


(Mb) 



- 1 . 


(25) 


(26) 


for every r, and that only r ef f (parameter a of gamma distribution) 
depends on r, while v eff (parameter b ) is constant for the whole profile. 
Then 


/ m , b m r[(\/b )—2 + m] w 
v ; ‘°p T[(l /b)—2] 

where 


(23) 


Let us consider the simplest case of a stratified cloud consisting 
of two layers. Let the top and bottom layers have droplet size 
distributions with common v e ff=b and different r e ff , respectively 
O] and a 2 (constant within each layer). The top layer has optical 
depth r,. 

9.2.1. Effective radius and variance 

In the 2-layer cloud case the w-moments (24) will have the form 


a m = | a m (T)w(7)dT, 
o 


(24) a m = a' 1 n jw( r) dr + a^ 1 j w(r)dr. 


(27) 
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Fig. 10. Summary of single layer simulations for r e ff=17.5 jjm and v e ff=0.1. Top: 
polarized reflectances vs. scattering angle for a variety of the layer optical depths. 
Middle: the residues of the reflectances from top left panel remaining after subtraction 
of the (scaled with the fitting factor v) reflectance corresponding to r=5“°° (solid 
lines), and linear fits to these residues (dashed lines). Bottom right: fit of u = (l — v) 
by an exponential function of t . 


Using the normalization condition (12), this expression can be 
written as 

0 ^= (a™—a™)W + a™, (28) 

where we introduced the notation of the top layer weight 


which can take values between 0 and 1. Note that the function w(r) 
can be different for different T\ (i.e., w = w(ti;t), thus, generally it 
is not a derivative of with respect to r-i). As follows from 

Eqs. (25) and (28), r e ff( top) can be expressed through W as 


Jtop) 


(a]-afjW + a\ 
+ aj ' 


(30) 


Note, that if a t = a 2 , the dependence on t , disappears, as expected. 
If the top layer is thin ( IV <k 1 ), we retrieve r^f° p) « a 2 , while r^f° p) « a, 
when it is thick (W«l). The weight W can be expressed using 
Eq. (30) in terms of the values of a,, a 2 , and r^° p) , and then substitut- 
ed back into Eqs. (28) and (26) to provide the relationship between 
the observed effective variance and the radius: 

v'“ p) =b-(1 + b)S, 6 2 , (31) 

where 


S f = 


.(top) 

etf 


-a, 


r 


(top) 

eff 


(32) 


for i = l,2 (see Appendix A for details of these calculations). Note, 
that since r^P p) is always in between a, and a 2 , S , and 8 2 have the op- 
posite signs. Thus, their product /ta, 6 2 < 0, and consequently v^fp p) > b. 
Note that v^r p) = b in the limit case of very thick or very thin top 
layer, when rif° p) is equal respectively to a , or a 2 (so that either 5] 
or d 2 is zero). As a function of r, (or, equivalently, of r£f° p) ) v£f° p) 
has its maximum value when r^rp 1 p) is equal to the harmonic mean 
of and a 2 : 


1 

,(“p) 


1(1 + 1 ) 
2 Vo, a 2 J 


(33) 


The maximal value of v£fp 1 p) is then 

max v'“ p) = b + (\ + b) ( 34 ) 

e " 4a, a 2 ' ’ 


9.2.2. Exponential weighting function 

The results presented so far are valid for weighting function w(r) 
of any shape. Let us now introduce a specific weighting function for a 
2-layer cloud following Eq. (14): 


w(r) = C m 


e _ftT 
e~ 


o-PiL o-ftCr-r,) 


r<r, 
r, <r<r„ 


(35) 


where C is the normalization constant, and s, and s 2 reflect the rela- 
tionship between droplet number concentration and optical depth 
(depending on particle size). Note that for gamma distribution 


5= [rr(Qext)a 2 (l-^)(l-2fa)] \ (36) 

and if we assume that <Q ext >«2 is the same for the two layers and 
so is the effective variance b, we can set 


Si 


1 a 1 

— 7T . and Sn — — rr . 

a] a 2 2 


(37) 


r, 

W = J w(r)dr, 
0 


(29) (Note, that the exact values of <Qext> can be incorporated 

by rescaling: a— *<Q ext > 1/2 a.) It is convenient to introduce the 
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notations u ^ and u 2 respectively for the top and bottom layer attenu- 
ation functions: 

U t =exp(-j3 1 r 1 ), (38) 

u 2 = exp[-p 2 (r max -Tj)]. (39) 


Using this notation 


C= £• + 


Ift [ft 

and therefore, 

1-Ui 




ftJ Ul 


W(r,) = 


where 




(40) 


(41) 


S = ?l-*(l-u 2 )=%L-%(l-u 2 ) 

Sl P 2 P 2 


.ft. 0? 


(42) 


Substitution of Eq. (41) into Eq. (30) yields the following expres- 
sion for the effective radius: 


jtop) 


= a, + 


(?) 

2 

£'( Q 2 — Q i) 

Q— T \ 

1 + 

Ks)V 

1 r l 


(43) 


Note that if T]=r max then f; = 0, and rif? p) = a,. If T] = 0, then 
exp( — = 1, and r^f° p) = a 2 . If the coefficients Sj are specified 
according to Eq, (37) and 0i =02 = j 3, then Eq. (43) takes the form 


r 


(top) 

eff 


= 0 , 


+ 1 - e -PTm„ 


(«2-ai): 


(44) 


which is especially simple in the case of semi-infinite medium: 

r eff P> = Q i + e ”' iri ( a 2-ai)- (45) 


The 2-layer distribution (35) can be naturally generalized to the 
case of a continuous profile, characterized by the dependences of 
s = tgff and the attenuation coefficient fi on the optical depth r: 


w(r) = Cs(r)exp 


-j(3(r)rdr 


(46) 


where the normalization constant C is determined from the condition 

( 12 ). 


9.2.3. Numerical results 

We conducted a series of numerical tests involving 2-layer clouds. 
In these tests simulated RSP measurements were made in the princi- 
pal plane for a SZA of 60". In all simulations the total optical depth of 
the 2-layer cloud was 5 (that is practically infinite for polarized re- 
flectance), while the top layer COD took the values of 0, 0.1, 0.2, 0.3, 
0.4, 0.5, 0.7, 1, 1.5, 2.5, 4, and 5. The layers had different r eff (7.5 and 
17.5 pm) and the same v e fp= 0.1. Two series of tests were conducted. 
In the first series (“Test A," Fig. 11 (left)) the layer with 17.5 pm drop- 
lets was on top, while in the second series (“Test B”, Fig. 11 (right)) 
the layer with smaller 7.5 pm droplets was on top. The RSP-type re- 
trievals of r£(f p) and v£fp p) are summarized in Table 1 and Fig. 12 (in 
the Figure, the nominal values of r e ff and v e ff were used instead of re- 
trievals at large 7j and at r, = 0). The maximum value of v eff allowed 
in the retrievals was 0.35. The top left plot in Fig. 12 shows that, as the 
top layer optical depth increases to Tj- 1-1.5, the retrieved value of 


r^ff p) reaches its limit value, which would be retrieved in the absence 
of the bottom layer. This convergence is apparently slower in Test A, 
where the top layer consists of larger particles. This means that 
even at the same optical depth the top layer with larger particles is 
more transparent than that with smaller ones, allowing for stronger 
influence of the bottom layer composition on the retrieved droplet 
size. This is consistent with the conjecture that small-angle multiple for- 
ward scattering events contribute to formation of the rainbow. Fig. 1 1 
shows that even when the observed polarized reflectance is a mixture 
of contributions from the 2 layers, and the effective size distribution 
n top is essentially bimodal (e.g., at r 1 =0.5), the Mie-theory-based LUT 
assuming a monomodal gamma distribution can still provide an ade- 
quate fit to the results, although a larger effective variance is retrieved. 
The bottom plots in Fig. 12 show that the width of the observed size dis- 
tribution first increases with the increase of the top layer optical depth 
(since now both top and bottom layers contribute to the observations), 
and then decreases as the influence of the bottom layer diminishes 
by Ti~ 1.5. The dotted lines in this plot show the values of v^ p) 
reconstructed from the corresponding r^r p) series using Eq. (31), which 
are in a reasonably good agreement with the direct numerical results. 
The theoretical maximum of vgr otect(top) , ; s 031 according to Eq. (34). It 
occurs (see Eq. (33)) at 4jr p) = 10.5 pm, that, corresponds to 7] = 0.34 
for Test A, and to 7j = 0.41 for Test B. We see that, as we already pointed 
out above, our direct retrievals of v e ff show a high bias. However, the 
positions of the maxima of the apparent effective variance in the direct 
retrievals are within 0.1 in optical depth of their theoretical estimates. 

Comparisons between analytical estimates and numerical results 
from Fig. 12 are presented in Fig. 13. The plots in this Figure compare 
the dependency on the top layer optical depth of the analytical esti- 
mates for the retrieved droplet effective radius r£fr p) , the retrieved 
effective variance and the top layer weights W(7f) with their numer- 
ical values obtained using Eq. (A.l). The bottom right plot in Fig. 13 
shows the relationship between numerical and analytical values of 
the large (17.5 pm) mode weight from the Tests A and B combined. 
While the values of the above described parameters retrieved from 
simulated data are in general agreement with the analytical predic- 
tions, certain systematic differences are clearly observed in Fig. 13. 
In particular, the numerical retrievals of r^° p) show slower than 
predicted change in the retrieved r£f° p) with increase of the top 
layer COD when the latter is small, and faster saturation when it is 
large. Also, the retrieved maximal values of v^fp p) are larger than 
predicted and occur at different top layer CODs. These discrepancies 
are consistent with the errors associated with fitting the results 
generated with a bi-modal size distribution using LUT computed as- 
suming that the distribution is mono-modal. Fig. 14 demonstrates 
these effects on fitting of both single-scattering polarized reflectances 
(left panels), and the size distribution shapes themselves (right 
panels). The retrievals here are presented as functions of the large 
(17.5 pm) mode statistical weight, thus, the plots of iif° p) and v£ff p) 
are similar to those in Fig. 13 corresponding to Test A. The bottom 
plots in Fig. 14 directly correspond to the bottom right plot in 
Fig. 13. All these plots show that least square fit algorithms tend to 
retrieve the parameters of the dominant mode in the bi-modal distri- 
bution rather than the averages over the whole distribution (that 
were computed analytically in these tests). 

10. Tests on 3D Monte Carlo simulations 

We used Monte Carlo simulations to study the influence of 3D ra- 
diative effects on RSP-type droplet size distribution retrievals. The ra- 
diative transfer model MYSTIC (Monte Carlo code for the phYSically 
correct Tracing of photons In Cloudy atmospheres; Mayer (2009)) 
was used for the simulation of RSP measurements for a 3D cloud 
field. MYSTIC is operated as one of several radiative transfer solvers 
of the libRadtran radiative transfer package (Mayer & Kylling, 2005). 
Originally, MYSTIC was developed as a forward tracing method for 
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Fig. 11. RSP-type retrievals from simulated polarized reflectance of 2-layer cloud with total COD = 5. Left: Test A (top layer r eff = 17.5 pm, bottom layer r eff = 7.5 pm). Right: Test B 
(top layer r e ff= 7.5 pm, bottom layer r e ff= 17.5 pm). In both tests v e tf= 0.1 for both layers. The top layer COD of is 0.1 (top), 0.5 (middle), and 1 (bottom). Red curves correspond to 
the actual multiple scattering simulations, while green curves show the best solution from the RSP retrieval algorithm. 


the calculation of irradiances and radiances in 3D plane-parallel at- 
mospheres. Recently, the model has been extended to include spher- 
ical geometry and a backward tracing mode (Emde & Mayer, 2007). In 
addition the model has been extended by Emde et al. (2010) to in- 
clude polarized radiation from scattering by randomly oriented parti- 
cles, i.e. clouds, aerosols, and molecules. MYSTIC includes variance 
reduction methods (Buras & Mayer, 2011) which are required for ef- 
ficient unbiased radiative transfer simulations in cloudy atmospheres. 

This model was applied to a realistic cloud field obtained from a 
large-eddy simulation (LES) of shallow, maritime convection. The 
LES model (Ackerman et al., 2004) treats three-dimensional fluid dy- 
namics of the atmosphere and incorporates a bin microphysics model 


that resolves the size distributions of aerosol and cloud droplets in 
each grid cell. It also includes a two-stream radiative transfer model 
that treats the vertical transport of radiation in each model column. 
The LES dataset used in this study is based on idealizations of mea- 
surements obtained during the Rain in Cumulus over the Ocean pro- 
ject (RICO, vanZanten et al., 2011), which was chosen to maximize 
the likelihood of observing 3D cloud effects on retrieved droplet 
sizes. This simulation has 100 m horizontal and 40 m vertical resolu- 
tion. To simplify the comparison between Monte Carlo computations 
and ID RT model, the actual bin-by-bin LES size distributions in 
each horizontal layer were replaced by gamma distributions with 
the same r eff , as in the microphysical model, and the effective variance 
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Table 1 

Results of RSP-type retrievals from simulated polarized reflectances of a 2-layer cloud 
with various top layer optical depth (total COD of both layers is 5). Test A: top layer 
r e(f =17.5 pm, bottom layer r eff =7.5 pm; Test B: vice versa. In all cases v eff =0.1 for 
both layers. 


Top r 

Test A 


Test B 


Top: r e ff=17.5 pm 
Bot.: r e ff=7.5 pm 


Top: r e ff = 7.5 pm 
Bot.: r e ff=17.5 pm 


r(ff p) , pm 

v‘ir° p > 

r<{? p) , pm 

v(ff P) 

0 

7.45 

0.117 

17.8 

0.100 

0.1 

7.55 

0.160 

16.8 

0.162 

0.2 

8.70 

0.285 

15.1 

0.285 

0.3 

9.90 

0.350 

12.2 

0.350 

0.4 

11.3 

0.350 

10.6 

0.350 

0.5 

13.5 

0.350 

9.65 

0.322 

0.7 

15.2 

0.292 

8.25 

0.192 

1 

16.7 

0.167 

7.50 

0.137 

1.5 

17.5 

0.117 

7.45 

0.120 

2.5 

17.8 

0.102 

7.50 

0.120 

4 

17.8 

0.100 

7.45 

0.117 

5 

17.8 

0.100 

7.45 

0.117 


uniformly set to 0.1. The simulations of RSP measurements were 
made assuming that the instrument is flown in the solar principal 
plane (coinciding with the y-axis of the LES grid) at 2.4 km above 
ground. The solar zenith angle used in the model run was 40°, the sur- 
face was modeled as a Lambertian reflector with a 5% albedo and the 
measurement wavelength was 554 nm. The US standard atmosphere 
was used for the pressure and trace gas profiles. 



TOP LAYER OPTICAL DEPTH 



To evaluate the effect of the 3D nature of the radiative fields as op- 
posed to the 1 D vertical structure, we also simulated the corresponding 
reflected radiation using an independent pixel approximation (IPA) 
(Cahalan et al., 1994) applied to the vertical profiles of cloud droplet 
sizes and concentrations provided by the LES model. A total of 35 cloudy 
pixels (optical depth greater than 0.5) along two transect lines 
(corresponding to two RSP-like “trajectories”) were originally selected 
for this comparison. The top left panel in Fig. 1 5 shows an example of 
the geometric shape of one of the cumulus clouds used. This shape is 
the cross-section along the transect line of points with droplet concen- 
trations greater than 20 cm -3 . The horizontal position of a “pixel” used 
for the comparison between ID and 3D retrievals is depicted by the 
dashed line. The RSP reflectance simulated with the 3D RT model output 
was aggregated to the point where this dashed line crosses the cloud 
top boundary. The corresponding ID simulations were made for a 
plane-parallel cloud with vertical profiles identical to those of the LES 
dataset along the dashed line. The profile of layer’s r eff is shown in the 
left middle panel of Fig. 15, while the corresponding 1 D and 3D simulat- 
ed polarized reflectance is presented in the top right panel together 
with the results of RSP-type droplet size retrievals (the corresponding 
fits are depicted by dotted curves). The values of size distribution pa- 
rameters retrieved for this pixel in the ID and 3D cases are close: rifp p) 
of respectively 13.2 and 13.5 pm, and v^f P) of respectively, 0.11 and 
0.12. These retrievals are consistent with r£ff p) = 13.75 pm computed 
directly from the LES model using a weighting function of the form 
given by Eq. (46). 

It is seen from the top right panel of Fig. 1 5 that the magnitude of 
3D polarized reflectance is substantially larger than that of its ID 




Fig. 12. Results of RSP retrievals made from simulated polarized reflectance of a 2-layer cloud. The data for Tests A and B are shown respectively in red and blue. Top left: depen- 
dence of the retrieved droplet effective radius rft?|i e ir on the top layer optical depth. The values of r e tf for each single layer are shown by dashed lines. Top right: The top layer 
weights W(r,) (from Eq. (A.l )). Bottom: dependence of the retrieved effective variance on the top layer optical depth. Dashed curves depict the values of v(Sf p) computed from 
r()f P) according to Eq. (31). 
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Fig. 13. Analytical estimations of the RSP retrievals in 2-layer cloud case (thick curves). The numerical results from Fig. 12 are shown by thin lines with diamonds depicting the 
actual data points. The data for Tests A and Bare shown respectively in red and blue. Top left: Dependence of the retrieved droplet effective radius ri“ p) on the top layer optical 
depth computed according to Eq. (30) with W from Eq. (41) and attenuation coefficients from Section 8. Top right: Retrieved effective variance computed according to Eq. (31) 
as function of the top layer optical depth. Bottom left: The top layer weights W(r,) analytical values are from Eq. (41) , the numerical values are obtained using Eq. (A.l). Bottom 
right: relationship between numerical and analytical values of the large (17.5 pm) mode weight from Tests A and B. Here the analytical weights are sampled at the data points of 
numerical simulations, the green curve represents a 5-degree polynomial fit to the numerical data. 


analog. This enhancement is due to the contribution of light entering 
the cloud from its sides in the 3D case. However, this enhancement 
has an effect only on the pixels on bright sides of clouds, which are 
exposed to the direct sunlight. For the pixels on the shadowed side 
of the clouds the opposite takes place: 3D simulations produce a 
smaller reflectance than the corresponding plane-parallel models. 
Since the polarized reflectance on the shadowed side of the cloud 
is often too small to exhibit a rainbow structure, we restrict further 
ID - 3D comparisons only to bright pixel cases. The middle right 
panel of Fig. 15 shows the histogram (for the set of all 35 cloudy 
pixels) of the ratio between 3D and ID polarized reflectances (taken 
at their minimum around 141"), where the bright pixel cases are 
separated from those in shadow by a dashed line. We see that for 
the most part bright pixel 3D polarized reflectance exceeds its ID 
counterpart by a factor of between 1.25 and 2. The bottom plots pres- 
ent comparisons between effective radius and variance retrievals 
from the ID and 3D radiative fields in 24 bright pixel cases. While 
no significant bias is seen between ID and 3D values of r£fp p) 
(0.2 pm mean difference and 0.73 pm standard deviation of the dif- 
ference), the corresponding v^f° p) values are systematically larger in 
the 3D case than in the ID case. The 3D values range within 0.13 ± 
0.02, while in the ID case v^fp p) is close to 0.1, which is the uniformly 
assumed single layer's v eff in the forward computations. This differ- 
ence can be explained by the larger penetration depth of light into 
the cloud in 3D cases, which results from the previously noted illumi- 
nation of the cloud from its sides. This effect causes a larger part of the 


non-uniform droplet size profile (cf. middle left panel) to contribute 
to the cloud-top retrievals, with a consequently wider droplet size 
distribution. 

11. Conclusions 

The presented sensitivity study of the accuracy of the cloud drop- 
let size retrievals from the RSP measurements demonstrated that this 
technique is a potentially valuable tool for climate research. 

Our algorithm retrieves the effective radius and variance of the cloud 
droplet size distribution, which is assumed to have a mono-modal 
gamma distribution shape. The tests were made on simulated datasets 
generated using both a plane-parallel (ID) and a Monte-Carlo (3D) RT 
model. We demonstrated that in the case of “measurements” made in 
the solar principle plane the accuracy of r cff retrievals is better on average 
than 0.15 pm (with the worst recorded case of 7.1 pm retrieval instead of 
7.5 pm at v eff = 0.2, that is a 5% error). Unlike r eff , the retrieved v eff shows 
systematic biases (see Fig. 4) ranging from 6% to 27% and generally 
decreasing with r e ir. The largest absolute differences occur at large v eff . 
This discrepancy can be explained by the “smoothing” effect on the 
reflectance of multiple scattering, which is interpreted by the retrieval 
algorithm as a larger v eff . The droplet size retrievals in an arbitrary 
viewing geometry require the Stokes vector to be transformed from the 
measurement coordinate frame into that of the scattering plane. In this 
case U is not identically zero and it is important to check that the param- 
eter Q, resulting from this reference frame rotation contains all the single 
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LARGE MODE WEIGHT - THEORY LARGE MODE WEIGHT - THEORY 


Fig. 14. Left: The effects of fitting of single-scattering polarized reflectance computed using Mie theory for bi-modal size distribution (with the same modes, as in Tests A and 
B) using LUT generated assuming mono-modal distribution model. Top: The retrieved effective radius (red) vs. the large (17.5 pm) mode weight in the bi-modal distribution. 
Black curve represents the analytical values from Eq. (30). Middle: Same as top but for effective variance. Bottom: large mode weight derived from the fitted r e fr (top plot) 
using Eq. (A.1) as function of the actual weight assumed in simulations. Right: same as left but for least square fit of size distributions themselves (rather than of derived optical 
parameters). 


scattering information on which our retrieval is based. We find that the 
difference between retrievals from rotated datasets and their principle 
plane counterparts is negligible: as little as 0.1-0.2 pm in radius and a 
few percentage points in variance. Our tests showed that uncertainties 
in scattering angle shift can be effectively resolved by considering the 
shift as a free parameter in the fitting procedure. This feature is 
implemented in our algorithm by introducing the shift parameter 6 in 
Eq. (3). We also determined that the shape of the rainbow in polarized 
reflectance is unchanged by overlaying aerosol scattering, although the 
signal is attenuated according to Beer's law. 

We demonstrated that the polarized rainbow amplitude includes 
contributions from both single and forward-directed multiple scattering. 


While the latter does not affect the droplet size retrievals, an assessment 
of the relative magnitudes of the various contributions is instrumental in 
understanding from where in a stratified cloud the remotely sensed radi- 
ative properties are generated. Our tests performed on a family of 2-layer 
clouds with widely different effective radii demonstrated that the param- 
eters retrieved from the simulated data are generally in agreement 
with the analytical predictions and show that for a typical airmass the 
retrieved size distribution parameters are generated within an optical 
depth ofl from cloud top (-50 m). There are however certain systematic 
differences, which indicate that the least square fit algorithm used tends 
to output the parameters of the dominant mode rather than an average 
over the whole bi-modal distribution. This is particularly relevant to 
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Fig. 15. Results of inter-comparison between droplet size retrievals from ID and 3D simulated radiation fields. Top left: shape of a cumulus cloud vertical cross-section from LES 
output. Top right: polarized reflectances and RSP-type droplet size retrievals corresponding to the cloud section depicted by the dashed line in the top left plot. Middle left: effective 
radius profile corresponding to this dashed line (from microphysical model). Middle right: frequency histogram of the ratio between 3D and ID polarized reflectances for 35 cloudy 
pixels. The dashed line separates parts of clouds illuminated by direct sunlight (bright cases) from those in shadow. Bottom left: comparison between effective radius values derived 
from ID and 3D fields (bright cases only). Bottom right: same as bottom left but for the effective variance. 


understanding under what conditions small evaporating drops near 
cloud top would be detected using polarization measurements. This 
problem may be potentially resolved by the Rainbow Fourier Transform 
(RFT), our new retrieval technique (Alexandrov et al., 2012) which al- 
lows for accurate retrieval of the size distribution shape of an arbitrary 
form (including multi-modal) with very good accuracy. 


Our final tests of the retrieval algorithm were applied to ID 1PA 
and 3D simulations of the same cloud field consisting of shallow con- 
vective clouds representative of those observed during RICO. The fact 
that 3D cloud effects can increase, or decrease the magnitude of the 
rainbow structure means that using a tabulation of exact ID RT 
results to reduce the bias in effective variance retrievals caused by 
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multiple scattering effects is unlikely to be effective. This is because 
multiple scattering effects are inherently sensitive to 3D effects. 
Nonetheless we find that the effective radius retrievals from the ID 
and 3D simulations agree extremely well while the effective variance 
retrieved from the 3D simulations is higher than that from the ID 
simulations. The reason for this is that for the simulated viewing ge- 
ometry (SZA 40 - in the principal plane) the rainbow is centered 
near nadir and side illumination of the clouds means that the rainbow 
is generated over a much larger depth into the clouds in the 3D sim- 
ulations than in the ID simulations. Vertical variations in effective 
radius are then interpreted as an increased effective variance. This re- 
sult is not “wrong" but emphasizes the importance of understanding 
the vertical weighting associated with any particular droplet size 
retrieval method (Platnick, 2000). Our plans for the future include 
testing our algorithm on a simulated cloud fields that include thin 
ice clouds above water clouds (this situation is common for measure- 
ments from space or high-altitude aircraft). 

In summary we find that the use of polarized reflectance observa- 
tions of the rainbow to retrieve the effective radius is extremely 
robust against both aerosol and 3D effects. The retrieved effective var- 
iance tends to be biased high in ID tests as a result of smoothing/ 
reduction in contrast of the rainbow structure caused by multiple 
scattering. The increased vertical depth over which the rainbow is 
generated in 3D clouds also tends to increase the retrieved effective 
variance. These biases are however relatively small and well under- 
stood and the retrieved effective variance is still a useful measure of 
droplet size dispersion near cloud top. 
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Appendix A. Relationship between effective radius and variance in 
2-layer cloud 


Eq. (30) can be used to express the top layer weight W through 
the observed effective radius and those of the two layers: 


W = 


hn 


fi 2 -hi ' 

where 


(A.1) 


fi i = ai i -afr7»= a f(a,.-r<r) 


(A.2) 


for i= 1,2. In these notations the w-moments (28) take the form: 

ho m 


h 2 -h, 
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In particular, 
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a 4 = J~rj^( a 2 — a l) [( a l + a 2) r eff PI ~ a l a 2 ] 

= ° 2 [( a l + a 2) r eff P — a l a 2 ] (A.6) 


When the expressions (A.4)-(A.6) are substituted into Eq. (26) for 
v<f? p) , the common factors a 2 = a|a|(a 2 — cii)/(h 2 — hi) cancel from the 
ratio and we obtain 


„(t°p) 


= (l + b) 


(Qi +a2)i'S P) -°i a 2 

( r eff P) ) 2 


- 1 , 


(A.7) 


that can also be written in the form of Eq. (31 ). 
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